#include <assert.h>
#include <stdio.h>

#include "../src/mstep.c"

static void
fill (double*           x,
      double          (*func) (void),
      unsigned int      n)
{
  for (unsigned int i = 0; i < n; ++i)
    {
      x[i] = func ();
    }
}

static void
fill_r (Rating*         x,
        unsigned int    n_workers,
        unsigned int    n_labels,
        unsigned int    n_images,
        unsigned int    n)
{
  for (unsigned int i = 0; i < n; ++i)
    {
      x[i].worker = (lrand48 () >> 16) % n_workers;
      x[i].image = (lrand48 () >> 16) % n_images;
      x[i].rating = (lrand48 () >> 16) % n_labels;
    }
}

static void
test_to_from_identity_once (void)
{
  /* (alpha, logbeta, gamma) = 
     from_x (alpha, logbeta, gamma, to_x (alpha, logbeta, gamma)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  double alpha[n_workers * n_labels * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double newalpha[n_workers * n_labels * n_labels];
  double newlogbeta[n_images];
  double newgamma[n_labels * n_labels];

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  fill (newalpha, drand48, n_workers * n_labels * n_labels);
  fill (logbeta, drand48, n_images);
  fill (newlogbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  fill (newgamma, drand48, n_labels * n_labels);

  size_t n_dims = 
    n_workers * n_labels * n_labels + n_images + n_labels * n_labels;
  gsl_vector* x = gsl_vector_alloc (n_dims);

  to_x (alpha, n_workers, n_labels, logbeta, n_images, gamma, x);
  from_x (newalpha, n_workers, n_labels, newlogbeta, n_images, newgamma, x);

  gsl_vector_free (x);

  assert (memcmp (alpha,
                  newalpha,
                  n_workers * n_labels * n_labels * sizeof (alpha[0])) == 0);

  assert (memcmp (logbeta, newlogbeta, n_images * sizeof (logbeta[0])) == 0);

  assert (memcmp (gamma, newgamma, n_labels * n_labels * sizeof (gamma[0])) == 0);
}

static void
test_to_from_inverse (void)
{
  for (unsigned int i = 0; i < 100; ++i)
    {
      test_to_from_identity_once ();
    }
}

static void
test_qfunc_f_once (void)
{
  /* qfunc (alpha, logbeta) = -qfunc_f (to_x (alpha, logbeta)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_ratings = ((lrand48 () >> 16) % 5) + 1;

  double alpha[n_workers * n_labels * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  Rating ratings[n_ratings];
  ParameterizedNormalDistribution prioralpha = parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution priorlogbeta = normal_distribution (drand48 (), 1.0 + drand48 ());
  NormalDistribution priorgamma = normal_distribution (drand48 (), 1.0 + drand48 ());

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  fill (logbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);
  fill_r (ratings, n_workers, n_labels, n_images, n_ratings);

  QfuncParams qfunc_params = { .n_workers = n_workers,
                               .n_labels = n_labels,
                               .n_images = n_images,
                               .pz = pz,
                               .priorz = priorz,
                               .ratings = ratings,
                               .n_ratings = n_ratings,
                               .prioralpha = &prioralpha.base,
                               .priorlogbeta = &priorlogbeta.base,
                               .priorgamma = &priorgamma.base,
                               .symmetric = 0,
                               .hyperprior = 1 };

  size_t n_dims = 
    n_workers * n_labels * n_labels + n_images + n_labels * n_labels;
  gsl_vector* x = gsl_vector_alloc (n_dims);

  to_x (alpha, n_workers, n_labels, logbeta, n_images, gamma, x);

  double q = qfunc (alpha, n_workers, n_labels, logbeta, n_images,
                    gamma, pz, priorz, ratings, n_ratings, 
                    &prioralpha.base, &priorlogbeta.base, &priorgamma.base,
                    0);

  double qf = qfunc_f (x, &qfunc_params);

  assert (fabs (q + qf) <= 1e-8 * (1 + fabs (q) + fabs (qf)) ||
          (fprintf (stderr, "%g ?= %g\n", q, qf), 0));

  gsl_vector_free (x);
}

static void
test_qfunc_f (void)
{
  for (unsigned int i = 0; i < 100; ++i)
    {
      test_qfunc_f_once ();
    }
}

static void
test_qfunc_df_once (void)
{
  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_ratings = ((lrand48 () >> 16) % 5) + 1;

  /* dqfunc (alpha, logbeta) = -qfunc_df (to_x (alpha, logbeta)) */

  double alpha[n_workers * n_labels * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double dalpha[n_workers * n_labels * n_labels];
  double dlogbeta[n_images];
  double dgamma[n_labels * n_labels];
  double newdalpha[n_workers * n_labels * n_labels];
  double newdlogbeta[n_images];
  double newdgamma[n_labels * n_labels];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  Rating ratings[n_ratings];
  ParameterizedNormalDistribution prioralpha = parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution priorlogbeta = normal_distribution (drand48 (), 1.0 + drand48 ());
  NormalDistribution priorgamma = normal_distribution (drand48 (), 1.0 + drand48 ());

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  fill (logbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  fill (dalpha, drand48, n_workers * n_labels * n_labels);
  fill (dlogbeta, drand48, n_images);
  fill (dgamma, drand48, n_labels * n_labels);
  fill (newdalpha, drand48, n_workers * n_labels * n_labels);
  fill (newdlogbeta, drand48, n_images);
  fill (newdgamma, drand48, n_labels * n_labels);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);
  fill_r (ratings, n_workers, n_labels, n_images, n_ratings);

  QfuncParams qfunc_params = { .n_workers = n_workers,
                               .n_labels = n_labels,
                               .n_images = n_images,
                               .pz = pz,
                               .priorz = priorz,
                               .ratings = ratings,
                               .n_ratings = n_ratings,
                               .prioralpha = &prioralpha.base,
                               .priorlogbeta = &priorlogbeta.base,
                               .priorgamma = &priorgamma.base,
                               .symmetric = 0,
                               .hyperprior = 1 };

  size_t n_dims = 
    n_workers * n_labels * n_labels + n_images + n_labels * n_labels;
  gsl_vector* x = gsl_vector_alloc (n_dims);
  gsl_vector* g = gsl_vector_alloc (n_dims);

  to_x (alpha, n_workers, n_labels, logbeta, n_images, gamma, x);

  qdqfunc (dalpha, dlogbeta, dgamma,
           alpha, n_workers, n_labels, logbeta, n_images, gamma,
           pz, priorz, ratings, n_ratings, 
           &prioralpha.base, &priorlogbeta.base, &priorgamma.base, 0);

  qfunc_df (x, &qfunc_params, g);

  from_x (newdalpha, n_workers, n_labels, newdlogbeta, n_images, newdgamma, g);

  for (unsigned int i = 0; i < n_workers * n_labels * n_labels; ++i)
    {
      assert (   fabs (newdalpha[i] + dalpha[i]) 
              <= 1e-8 * (1 + fabs (newdalpha[i]) + fabs (dalpha[i])) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, newdalpha[i], dalpha[i]),
               0));
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (   fabs (newdlogbeta[i] + dlogbeta[i]) 
              <= 1e-8 * (1 + fabs (newdlogbeta[i]) + fabs (dlogbeta[i])) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, newdlogbeta[i], dlogbeta[i]),
               0));
    }

  for (unsigned int i = 0; i < n_labels * n_labels; ++i)
    {
      assert (   fabs (newdgamma[i] + dgamma[i]) 
              <= 1e-8 * (1 + fabs (newdgamma[i]) + fabs (dgamma[i])) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, newdgamma[i], dgamma[i]),
               0));
    }

  gsl_vector_free (g);
  gsl_vector_free (x);
}

static void
test_qfunc_df (void)
{
  for (unsigned int i = 0; i < 100; ++i)
    {
      test_qfunc_df_once ();
    }
}

static void
test_mstep_no_data_once (void)
{
  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  /* without data mstep should converge on prior maximum */

  double alpha[n_workers * n_labels * n_labels];
  double diagonal[n_workers * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double gdiagonal[n_labels];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  ParameterizedNormalDistribution prioralpha = parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution priorlogbeta = normal_distribution (drand48 (), 1.0 + drand48 ());
  NormalDistribution priorgamma = normal_distribution (drand48 (), 1.0 + drand48 ());

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  fill (logbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);

  for (unsigned int j = 0; j < n_workers; ++j)
    {
      for (unsigned int k = 0; k < n_labels; ++k)
        {
          diagonal[j * n_labels + k] = 
            alpha[j * n_labels * n_labels + k * n_labels + k];
          gdiagonal[k] = gamma[k * n_labels + k];
        }
    }

  mstep (alpha, n_workers, n_labels, 
         logbeta, n_images, gamma,
         pz, priorz, 
         NULL, 0, 
         &prioralpha.base, &priorlogbeta.base, &priorgamma.base,
         0, 1,
         50, 1e-8);

  for (unsigned int i = 0, x = 0; i < n_workers * n_labels * n_labels; ++i)
    {
      if (i % n_labels == (i / n_labels) % n_labels)
        {
          assert (   fabs (alpha[i] - diagonal[x]) <= 
                     1e-8 * (1 + fabs (alpha[i]) + fabs (diagonal[x])) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, alpha[i], diagonal[x]),
                   0));

          ++x;
        }
      else
        {
          assert (   fabs (alpha[i] - priorgamma.mean) <= 
                     1e-8 * (1 + fabs (alpha[i]) + fabs (priorgamma.mean)) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, alpha[i], priorgamma.mean),
                   0));
        }
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (   fabs (logbeta[i] - priorlogbeta.mean) <= 
                 1e-8 * (1 + fabs (logbeta[i]) + fabs (priorlogbeta.mean)) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, logbeta[i], priorlogbeta.mean),
               0));
    }

  for (unsigned int i = 0, x = 0; i < n_labels * n_labels; ++i)
    {
      if (i % n_labels == (i / n_labels) % n_labels)
        {
          assert (   fabs (gamma[i] - gdiagonal[x]) <= 
                     1e-8 * (1 + fabs (gamma[i]) + fabs (gdiagonal[x])) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], gdiagonal[x]),
                   0));

          ++x;
        }
      else
        {
          assert (   fabs (gamma[i] - priorgamma.mean) <= 
                     1e-8 * (1 + fabs (gamma[i]) + fabs (priorgamma.mean)) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], priorgamma.mean),
                   0));
        }
    }
}

static void
test_mstep_no_data (void)
{
  for (unsigned int i = 0; i < 20; ++i)
    {
      test_mstep_no_data_once ();
    }
}

static void
test_mstep_symmetric_no_data_once (void)
{
  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  /* without data mstep should converge on prior maximum */

  double alpha[n_workers * n_labels * n_labels];
  double alphacopy[n_workers * n_labels * n_labels];
  double logbeta[n_images];
  double gamma[n_labels * n_labels];
  double gammacopy[n_labels * n_labels];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  ParameterizedNormalDistribution prioralpha = parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution priorlogbeta = normal_distribution (drand48 (), 1.0 + drand48 ());
  NormalDistribution priorgamma = normal_distribution (drand48 (), 1.0 + drand48 ());

  fill (alpha, drand48, n_workers * n_labels * n_labels);
  memcpy (alphacopy, alpha, n_workers * n_labels * n_labels * sizeof (alpha[0]));
  fill (logbeta, drand48, n_images);
  fill (gamma, drand48, n_labels * n_labels);
  memcpy (gammacopy, gamma, n_labels * n_labels * sizeof (gamma[0]));
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);

  mstep (alpha, n_workers, n_labels, 
         logbeta, n_images, gamma,
         pz, priorz, 
         NULL, 0, 
         &prioralpha.base, &priorlogbeta.base, &priorgamma.base, 
         1, 1,
         50, 1e-8);

  for (unsigned int i = 0; i < n_workers * n_labels * n_labels; ++i)
    {
      if (i % n_labels <= (i / n_labels) % n_labels)
        {
          assert (   fabs (alpha[i] - alphacopy[i]) <= 
                     1e-8 * (1 + fabs (alpha[i]) + fabs (alphacopy[i])) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, alpha[i], alphacopy[i]),
                   0));
        }
      else
        {
          assert (   fabs (alpha[i] - priorgamma.mean) <= 
                     1e-8 * (1 + fabs (alpha[i]) + fabs (priorgamma.mean)) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, alpha[i], priorgamma.mean),
                   0));
        }
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (   fabs (logbeta[i] - priorlogbeta.mean) <= 
                 1e-8 * (1 + fabs (logbeta[i]) + fabs (priorlogbeta.mean)) ||
              (fprintf (stderr, "%u: %g ?= %g\n", i, logbeta[i], priorlogbeta.mean),
               0));
    }

  for (unsigned int i = 0; i < n_labels * n_labels; ++i)
    {
      if (i % n_labels <= (i / n_labels) % n_labels)
        {
          assert (   fabs (gamma[i] - gammacopy[i]) <= 
                     1e-8 * (1 + fabs (gamma[i]) + fabs (gammacopy[i])) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], gammacopy[i]),
                   0));
        }
      else
        {
          assert (   fabs (gamma[i] - priorgamma.mean) <= 
                     1e-8 * (1 + fabs (gamma[i]) + fabs (priorgamma.mean)) ||
                  (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], priorgamma.mean),
                   0));
        }
    }
}

static void
test_mstep_symmetric_no_data (void)
{
  for (unsigned int i = 0; i < 20; ++i)
    {
      test_mstep_symmetric_no_data_once ();
    }
}

typedef struct _MstepTestDatum MstepTestDatum ;
struct _MstepTestDatum 
{ 
  struct
    {
      double    alpha[4 * 5 * 5];
      double    logbeta[3];
      double    gamma[5 * 5];
      double    pz[3 * 5];
      double    priorz[5];
      Rating    ratings[8];
      double    alphastddev;
      double    logbetamean;
      double    logbetastddev;
      double    gammamean;
      double    gammastddev;
    }                   inputs;

  struct
    {
      double    newq;
      struct
        {
          double    alpha[4 * 5 * 5];
          double    logbeta[3];
          double    gamma[5 * 5];
        }       new;
    }                   desired_output;
};

static void
test_mstep (void)
{
  MstepTestDatum data[] = 
    {{{{0, 0.290552, 0.211452, 0.2384, 0.0482029, 0.603504, 0, 0.0200152,  0.0490884, 0.581799, 0.129669, 0.794694, 0, 0.039053, 0.693316,  0.890465, 0.503541, 0.0560458, 0, 0.975266, 0.516066, 0.687009,  0.0910183, 0.88197, 0, 0, 0.77709, 0.910482, 0.342099, 0.425496,  0.486538, 0, 0.69903, 0.103699, 0.377293, 0.883034, 0.679015, 0,  0.0546107, 0.795494, 0.753365, 0.884321, 0.0155576, 0, 0.102178,  0.8629, 0.38078, 0.959512, 0.126913, 0, 0, 0.346834, 0.69377,  0.868494, 0.244943, 0.569744, 0, 0.783288, 0.526395, 0.819447,  0.0832067, 0.0842582, 0, 0.422696, 0.442154, 0.200173, 0.405243,  0.368085, 0, 0.64666, 0.446808, 0.520922, 0.352527, 0.544481, 0,  0, 0.583908, 0.140143, 0.393016, 0.417568, 0.237074, 0, 0.446372,  0.524522, 0.172626, 0.66733, 0.663084, 0, 0.998127, 0.353179,  0.584123, 0.578826, 0.575432, 0, 0.911025, 0.383951, 0.173582,  0.207346, 0.264365, 0}, {0.937143, 0.65266, 0.854819}, {1.,  0.719884, 0.353234, 0.512518, 0.461803, 0.302315, 1., 0.11616,  0.0661453, 0.937281, 0.12969, 0.44883, 1., 0.403061, 0.939154,  0.776511, 0.864707, 0.824236, 1., 0.363722, 0.865486, 0.480756,  0.650653, 0.156376, 1.}, {0.601121, 0.543614, 0.997993, 0.301557,  0.881238, 0.19038, 0.485476, 0.839754, 0.578922, 0.0742193,  0.41933, 0.902473, 0.449233, 0.625389, 0.0162689}, {0.963319,  0.672722, 0.760682, 0.192033,  0.599597}, {{1, 1, 2}, {0, 2, 0}, {0, 1, 1}, {3, 1, 0}, {3, 1,  2}, {3, 0, 1}, {0, 1, 3}, {3, 0, 0}}, 0.807236, 0.279926,  0.54138, 0.443221,  0.206115}, {-121.419, {{0., 0.2771, 0.549331, 0.294221, 0.576189,  0.212499, 0., 0.626759, 0.163621, 0.671747, 0.364645, -0.080317,  0., -0.0729796, 0.642284, 0.330252, 0.113418, 0.617561, 0.,  0.666759, 0.437439, 0.320358, 0.476308, 0.354665, 0., 0.,  0.463271, 0.16397, 0.486471, 0.501684, 0.480127, 0., -0.105163,  0.510986, 0.536305, 0.586703, 0.582153, 0., 0.602279, 0.631049,  0.490407, 0.500852, -0.165369, 0., 0.539387, 0.432427, 0.42918,  0.325783, 0.466404, 0., 0., 0.423377, 0.432483, 0.449207,  0.466049, 0.405434, 0., 0.407434, 0.442769, 0.472981, 0.416613,  0.409772, 0., 0.4398, 0.481746, 0.410839, 0.423689, 0.396634, 0.,  0.470488, 0.414311, 0.41089, 0.443166, 0.45001, 0., 0.,  0.225389, 0.419455, 0.658757, 0.670435, -0.0559206, 0., 0.151801,  0.64678, 0.667377, -0.109631, 0.214437, 0., 0.737639,  0.762836, -0.0848313, 0.35722, 0.0231363, 0.,  0.623557, -0.0703569, -0.012771, 0.526568, 0.633102,  0.}, {0.63372, 1.07827, 0.26912}, {1., 0.423377, 0.432483,  0.449207, 0.466049, 0.405434, 1., 0.407434, 0.442769, 0.472981,  0.416613, 0.409772, 1., 0.4398, 0.481746, 0.410839, 0.423689,  0.396634, 1., 0.470488, 0.414311, 0.41089, 0.443166, 0.45001,  1.}}}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution prior_alpha = 
        parameterized_normal_distribution (data[i].inputs.alphastddev);
      NormalDistribution prior_logbeta = 
        normal_distribution (data[i].inputs.logbetamean, 
                             data[i].inputs.logbetastddev);
      NormalDistribution prior_gamma = 
        normal_distribution (data[i].inputs.gammamean, 
                             data[i].inputs.gammastddev);
      double diagonal[4 * 5];
      double gdiagonal[5];

      for (unsigned int j = 0; j < 4; ++j)
        {
          for (unsigned int k = 0; k < 5; ++k)
            {
              diagonal[j * 5 + k] = drand48 ();
              data[i].inputs.alpha[j * 5 * 5 + k * 5 + k] = diagonal[j * 5 + k];
            }
        }

      for (unsigned int k = 0; k < 5; ++k)
        {
          gdiagonal[k] = drand48 ();
          data[i].inputs.gamma[k * 5 + k] = gdiagonal[k];
        }

      double q = mstep (data[i].inputs.alpha,
                        4,
                        5,
                        data[i].inputs.logbeta,
                        3,
                        data[i].inputs.gamma,
                        data[i].inputs.pz,
                        data[i].inputs.priorz,
                        data[i].inputs.ratings,
                        8,
                        &prior_alpha.base,
                        &prior_logbeta.base,
                        &prior_gamma.base,
                        0,
                        1,
                        100,
                        1e-8);

      assert (fabs (q - data[i].desired_output.newq) <=
              1e-4 * (1 + fabs (q) + fabs (data[i].desired_output.newq)) ||
              ( fprintf (stderr, "%g ?= %g\n", q, data[i].desired_output.newq),
                0 ));
      
      for (unsigned int j = 0, y = 0; j < 4*5*5; ++j)
        {
          if (j % 5 == (j / 5) % 5)
            {
              assert (fabs (data[i].inputs.alpha[j] - diagonal[y]) <=
                      1e-8 * (1 + fabs (data[i].inputs.alpha[j]) + fabs (diagonal[y])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.alpha[j], diagonal[y]),
                        0 ));

              ++y;
            }
          else
            {
              assert (fabs (data[i].inputs.alpha[j] - data[i].desired_output.new.alpha[j]) <=
                      5e-3 * (1 + fabs (data[i].inputs.alpha[j]) + fabs (data[i].desired_output.new.alpha[j])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.alpha[j], data[i].desired_output.new.alpha[j]),
                        0 ));
            }
        }

      for (unsigned int j = 0; j < 3; ++j)
        {
          assert (fabs (data[i].inputs.logbeta[j] - data[i].desired_output.new.logbeta[j]) <=
                  5e-3 * (1 + fabs (data[i].inputs.logbeta[j]) + fabs (data[i].desired_output.new.logbeta[j])) ||
                  ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.logbeta[j], data[i].desired_output.new.logbeta[j]),
                    0 ));
        }

      for (unsigned int j = 0, y = 0; j < 5*5; ++j)
        {
          if (j % 5 == (j / 5) % 5)
            {
              assert (fabs (data[i].inputs.gamma[j] - gdiagonal[y]) <=
                      1e-8 * (1 + fabs (data[i].inputs.gamma[j]) + fabs (diagonal[y])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.gamma[j], diagonal[y]),
                        0 ));

              ++y;
            }
          else
            {
              assert (fabs (data[i].inputs.gamma[j] - data[i].desired_output.new.gamma[j]) <=
                      5e-3 * (1 + fabs (data[i].inputs.gamma[j]) + fabs (data[i].desired_output.new.gamma[j])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.gamma[j], data[i].desired_output.new.gamma[j]),
                        0 ));
            }
        }
    }
}

typedef struct _MstepNoHyperTestDatum MstepNoHyperTestDatum ;
struct _MstepNoHyperTestDatum 
{ 
  struct
    {
      double    alpha[4 * 5 * 5];
      double    logbeta[3];
      double    gamma[5 * 5];
      double    pz[3 * 5];
      double    priorz[5];
      Rating    ratings[8];
      double    alphastddev;
      double    logbetamean;
      double    logbetastddev;
      double    gammamean;
      double    gammastddev;
    }                   inputs;

  struct
    {
      double    newq;
      struct
        {
          double    alpha[4 * 5 * 5];
          double    logbeta[3];
        }       new;
    }                   desired_output;
};

static void
test_mstep_no_hyperprior (void)
{
  MstepNoHyperTestDatum data[] = 
    {{{{0, 0.892916, 0.757276, 0.889142, 0.239236, 0.541408, 0, 0.864213,  0.791336, 0.729487, 0.476928, 0.825879, 0, 0.268104, 0.575471,  0.770123, 0.988367, 0.891043, 0, 0.787295, 0.741929, 0.55213,  0.839032, 0.782372, 0, 0, 0.671842, 0.0395729, 0.255145, 0.454443,  0.778925, 0, 0.282297, 0.366003, 0.215207, 0.237517, 0.418084, 0,  0.574667, 0.48572, 0.760589, 0.592204, 0.306563, 0, 0.910249,  0.990466, 0.603837, 0.41552, 0.122954, 0, 0, 0.248537, 0.0517074,  0.576487, 0.340582, 0.576696, 0, 0.0121345, 0.321342, 0.886139,  0.79777, 0.729838, 0, 0.955339, 0.670932, 0.560253, 0.311754,  0.380672, 0, 0.185212, 0.799664, 0.71955, 0.0741091, 0.274963, 0,  0, 0.809198, 0.115712, 0.658589, 0.152009, 0.560661, 0, 0.0640051,  0.0821016, 0.811427, 0.983965, 0.0518706, 0, 0.76076, 0.925288,  0.186195, 0.322033, 0.805421, 0, 0.254355, 0.625941, 0.0102784,  0.424749, 0.0691433, 0}, {0.826277, 0.290729, 0.350639}, {1.,  0.79418, 0.0170788, 0.175016, 0.69205, 0.642171, 1., 0.456418,  0.111011, 0.609949, 0.830745, 0.472453, 1., 0.0591404, 0.849189,  0.905457, 0.286259, 0.737108, 1., 0.0437687, 0.651102, 0.660317,  0.726829, 0.61902, 1.}, {0.581958, 0.83404, 0.436101, 0.268381,  0.787778, 0.816961, 0.261085, 0.57633, 0.145607, 0.360543,  0.150074, 0.966381, 0.314862, 0.88809, 0.0909332}, {0.117192,  0.409405, 0.601831, 0.353825,  0.0734236}, {{2, 0, 2}, {3, 1, 0}, {0, 1, 2}, {3, 0, 3}, {0, 1,  3}, {1, 1, 3}, {0, 1, 0}, {3, 0, 0}}, 0.758303, 0.941514,  0.626996, 0.454404,  0.176345}, {-150.71, {{0., 0.869154, 0.0012881, 0.0249056,  0.791112, 0.304716, 0., 0.20442, 0.0565133, 0.66483, 0.25726,  0.632271, 0., 0.0691535, 0.908752, 0.652065, 0.398688, 0.512975,  0., 0.241102, 0.286272, 0.735741, 0.319739, 0.27275, 0., 0.,  0.811903, 0.183884, -0.355977, 0.716927, 0.659175, 0.,  0.487772, -0.135156, 0.628894, 0.844719, 0.517408, 0., -0.297987,  0.862313, 0.911143, 0.331315, 0.747313, 0., 0.13468, 0.685371,  0.693594, 0.753682, 0.128814, 0., 0., 0.797793, -0.321825,  0.244863, 0.698209, 0.657618, 0., -0.238224, 0.248919, 0.628046,  0.84213, 0.532307, 0., 0.284664, 0.85956, 0.908568, 0.348679,  0.121372, 0., 0.189892, 0.672586, 0.680869, -0.128695, 0.644052,  0., 0., 0.855614, 0.423485, 0.112779, 0.777552, -0.165207, 0.,  0.525287, -0.164837, 0.647821, -0.106275, 0.549583,  0., -0.0533938, 0.869112, 0.195409, 0.420816, 0.76285, 0.,  0.2952, -0.111623, 0.702962, 0.759555, -0.0505377, 0.}, {1.67833,  1.27253, 0.941514}}}}}; 

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution prior_alpha = 
        parameterized_normal_distribution (data[i].inputs.alphastddev);
      NormalDistribution prior_logbeta = 
        normal_distribution (data[i].inputs.logbetamean, 
                             data[i].inputs.logbetastddev);
      NormalDistribution prior_gamma = 
        normal_distribution (data[i].inputs.gammamean, 
                             data[i].inputs.gammastddev);
      double diagonal[4 * 5];
      double gammacopy[5 * 5];

      for (unsigned int j = 0; j < 4; ++j)
        {
          for (unsigned int k = 0; k < 5; ++k)
            {
              diagonal[j * 5 + k] = drand48 ();
              data[i].inputs.alpha[j * 5 * 5 + k * 5 + k] = diagonal[j * 5 + k];
            }
        }

      memcpy (gammacopy, data[i].inputs.gamma, 5 * 5 * sizeof (gammacopy[0]));

      double q = mstep (data[i].inputs.alpha,
                        4,
                        5,
                        data[i].inputs.logbeta,
                        3,
                        data[i].inputs.gamma,
                        data[i].inputs.pz,
                        data[i].inputs.priorz,
                        data[i].inputs.ratings,
                        8,
                        &prior_alpha.base,
                        &prior_logbeta.base,
                        &prior_gamma.base,
                        0,
                        0,
                        100,
                        1e-8);

      assert (fabs (q - data[i].desired_output.newq) <=
              1e-4 * (1 + fabs (q) + fabs (data[i].desired_output.newq)) ||
              ( fprintf (stderr, "%g ?= %g\n", q, data[i].desired_output.newq),
                0 ));
      
      for (unsigned int j = 0, y = 0; j < 4*5*5; ++j)
        {
          if (j % 5 == (j / 5) % 5)
            {
              assert (fabs (data[i].inputs.alpha[j] - diagonal[y]) <=
                      1e-8 * (1 + fabs (data[i].inputs.alpha[j]) + fabs (diagonal[y])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.alpha[j], diagonal[y]),
                        0 ));

              ++y;
            }
          else
            {
              assert (fabs (data[i].inputs.alpha[j] - data[i].desired_output.new.alpha[j]) <=
                      5e-3 * (1 + fabs (data[i].inputs.alpha[j]) + fabs (data[i].desired_output.new.alpha[j])) ||
                      ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.alpha[j], data[i].desired_output.new.alpha[j]),
                        0 ));
            }
        }

      for (unsigned int j = 0; j < 3; ++j)
        {
          assert (fabs (data[i].inputs.logbeta[j] - data[i].desired_output.new.logbeta[j]) <=
                  5e-3 * (1 + fabs (data[i].inputs.logbeta[j]) + fabs (data[i].desired_output.new.logbeta[j])) ||
                  ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.logbeta[j], data[i].desired_output.new.logbeta[j]),
                    0 ));
        }

      for (unsigned int j = 0; j < 5*5; ++j)
        {
          assert (fabs (data[i].inputs.gamma[j] - gammacopy[j]) <=
                  1e-8 * (1 + fabs (data[i].inputs.gamma[j]) + fabs (gammacopy[j])) ||
                  ( fprintf (stderr, "%u %u: %g ?= %g\n", i, j, data[i].inputs.gamma[j], gammacopy[j]),
                    0 ));
        }
    }
}

int 
main (void)
{
  srand48 (69);

  test_to_from_inverse ();
  test_qfunc_f ();
  test_qfunc_df ();
  test_mstep_no_data ();
  test_mstep_symmetric_no_data ();
  test_mstep ();
  test_mstep_no_hyperprior ();

  return 0;
}
